kinfitr Function Examples

Granville J Matheson

2017-01-04

Background

Reproducible research involves analysis and reporting of the data from raw data to final outcomes in one unified framework. Given the data, another researcher should be able to run it through the same analysis, and see each step of the analysis and check the robustness of the findings as they wish. This also means that the addition of extra data, or omission of problematic data, or changes to the way that data should be processed should be straightforward to implement without the need of error-prone manual entering and reporting of data. In this way, research is made more transparent, errors should be rarer and easier to diagnose, and the analysis itself is also made easier in the long run, with its requiring little to no extra work as data sets increase, or errors in the raw data are discovered.

R is the ideal option for performing reproducible research as its primary IDE, RStudio, comes prepackaged with R Markdown and knitr, which is specifically designed for reproducible research and generating reproducible reports. Furthermore, as opposed to MATLAB, it is free which means that the data analysis can be run by anyone and not only those who pay for the software. It is also a language designed for statistical analysis, meaning that it is ideally suited to the statistical analysis and presentation. Finally, there are powerful tools within R for performing the fitting procedures required for effective kinetic modelling.

Software available at present for TAC extraction and kinetic analysis include PMOD (which is expensive and which cannot be scripted as far as I know), MIAKAT (which is open-source, but still requires MATLAB which is not freeware), and various custom-built tools within groups, most often built in MATLAB. Statistical modelling and presentation are sometimes performed in MATLAB, but are usually performed using various statistical and graphical programs.

All the processes up until TAC extraction require large raw data files to be shared and to be reproducible, which are often complicated not only to store, but also to share as they contain sensitive patient information in the form of an image, which is theoretically difficult to fully anonymise and therefore ethically problematic. TAC extraction is also a process which can be performed in many different software packages as this does not differ between PET and fMRI, including, for example, NiPy in Python, which is also a free software package. Storage and sharing of TACs is more convenient as this data is in the form of a text document containing one-dimensional time series data. This therefore, in the opinion of the author, represents the ideal step of the image analysis pipeline to begin a reproducible analysis which can be shared in full. This can also be kept for documentation purposes in the interests of data transparency.

For these reasons, and to enhance my own knowledge of popular kinetic modelling procedures, I have written a kinetic modelling R package which can analyse and plot TAC data using several popular kinetic models, while easily accepting input data for each function in as raw a form as possible without requiring complex data structures. Furthermore, each process stores not only the parameters, but also the fit itself, the fitted values, all of the input data, and has straightforward tools for plotting the model to perform quality control.

# Loading packages
library(kinfitr)

library(knitr)
library(pander)
library(broom)

Reference Region Models

Data

times = tacdata$Time/60
reftac <- tacdata$RefCBL
roitac <- tacdata$FSLSSTR
weights <- tacdata$Weights

First all the data is loaded (unshown). This data includes:

This data is plotted below:

par(mfrow=c(3,1))

plot(times, reftac)
plot(times, roitac)
plot(times, weights)
Reference Region Model Data

Reference Region Model Data

SRTM

\[C(t) = R_{1}C_{R}(t)+(k_{2}-\frac{R_{1}k_{2}}{1+BP_{ND}})C_{R}(t)\otimes e^{-\frac{k_{2}}{1+BP_{ND}}t}\]

SRTM requires input of the times, reftac and roitac. Weights are optional.

srtmout <- srtm(times, reftac, roitac, weights)

Parameters

pandoc.table(srtmout$par, digits = 3)
R1 k2 bp
0.856 0.111 1.53

Plot

The output object can be used to visualise the fit:

plot_srtmfit(srtmout)
SRTM Fit Check

SRTM Fit Check

… and the residuals:

plot_residuals(srtmout)
SRTM Residuals

SRTM Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(srtmout))
  Length Class Mode
par 3 data.frame list
par.se 3 data.frame list
fit 6 nls list
weights 14 -none- numeric
tacs 4 data.frame list
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(srtmout$fit), digits=3)
sigma isConv finTol logLik AIC BIC deviance df.residual
5.43 TRUE 1.49e-08 -57.4 123 125 295 10

Non-Invasive Logan Plot

refLogan requires input of the times, reftac, roitac and k2prime. Weights are optional.

k2prime can be obtained from the srtm fit:

k2prime = srtmout$par$k2 / srtmout$par$R1
k2prime
## [1] 0.1292453

First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.

taclow = tacdata$FSLSFC    # Frontal Cortex 
tacmed = tacdata$fslLST    # Limbic Striatum
tachigh = tacdata$fslAST   # Associative Striatum

refLogan_tstar(times, reftac, taclow, tacmed, tachigh, k2prime = k2prime, 'demonstration')
## png 
##   2

refLogan tstar finder

Using this information, we see that we should probably use 9 included frames for the fit.

refloganout <- refLogan(times, reftac, roitac, k2prime, 
                        tstarIncludedFrames = 9, weights=weights)

Parameters

pandoc.table(refloganout$par)
bp
1.522

Plot

The output object can be used to visualise the fit:

plot_refLoganfit(refloganout)
RefLogan Fit Check

RefLogan Fit Check

… and the residuals:

plot_residuals(refloganout)
RefLogan Residuals

RefLogan Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(refloganout))
  Length Class Mode
par 1 data.frame list
fit 13 lm list
tacs 3 data.frame list
fitvals 2 data.frame list
weights 14 -none- numeric
k2prime 1 -none- numeric
tstarIncludedFrames 1 -none- numeric
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(refloganout$fit), digits = 3)
Table continues below
r.squared adj.r.squared sigma statistic p.value df logLik AIC
1 1 0.279 40528 1.97e-14 2 -0.742 7.48
BIC deviance df.residual
8.08 0.543 7

MRTM1

MRTM1 requires input of the times, reftac and roitac. Weights are optional.

First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.

mrtm1_tstar(times, reftac, taclow, tacmed, tachigh, 'demonstration')

Using this information (unshown), we see that we should probably use all 14 frames for the fit.

mrtm1out <- mrtm1(times, reftac, roitac, tstarIncludedFrames = 14, weights=weights)

Parameters

pandoc.table(mrtm1out$par, digits = 3)
bp k2prime
1.52 0.13

Plot

The output object can be used to visualise the fit:

plot_mrtm1fit(mrtm1out)
MRTM1 Fit Check

MRTM1 Fit Check

… and the residuals:

plot_residuals(mrtm1out)
MRTM1 Residuals

MRTM1 Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(mrtm1out))
  Length Class Mode
par 2 data.frame list
fit 13 lm list
tacs 3 data.frame list
fitvals 9 data.frame list
weights 14 -none- numeric
tstarIncludedFrames 1 -none- numeric
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(mrtm1out$fit), digits = 3)
Table continues below
r.squared adj.r.squared sigma statistic p.value df logLik AIC
1 1 5.1 16740 8.46e-19 3 -39.2 86.4
BIC deviance df.residual
88.7 260 10

MRTM2

MRTM2 requires input of the times, reftac, roitac and k2prime. Weights are optional.

k2prime can be obtained from the mrtm1 fit:

k2prime = mrtm1out$par$k2prime
k2prime
## [1] 0.1295072

First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.

mrtm2_tstar(times, reftac, taclow, tacmed, tachigh, k2prime=k2prime, 'demonstration')

Using this information (unshown), we see that we should probably use all 14 frames for the fit.

mrtm2out <- mrtm2(times, reftac, roitac, k2prime = k2prime, 
                  tstarIncludedFrames = 14, weights=weights)

Parameters

pandoc.table(mrtm2out$par)
bp
1.521

Plot

The output object can be used to visualise the fit:

plot_mrtm2fit(mrtm2out)
MRTM2 Fit Check

MRTM2 Fit Check

… and the residuals:

plot_residuals(mrtm2out)
MRTM2 Residuals

MRTM2 Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(mrtm2out))
  Length Class Mode
par 1 data.frame list
fit 13 lm list
tacs 3 data.frame list
fitvals 8 data.frame list
weights 14 -none- numeric
k2prime 1 -none- numeric
tstarIncludedFrames 1 -none- numeric
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(mrtm2out$fit), digits = 3)
Table continues below
r.squared adj.r.squared sigma statistic p.value df logLik AIC
1 1 4.86 27622 4.41e-21 2 -39.2 84.4
BIC deviance df.residual
86.1 260 11

Non-Invasive Multilinear Logan Plot

\[\int _{0} ^{t} C_{T} (\tau) d\tau = DVR\Bigg(\int _{0} ^{t} C_{T}' (\tau) d\tau+\frac{C_{T}'(t)}{k2'}\Bigg)-bC_{T}(t)\]

refmlLogan is the kinetic model used in WAPI, which is detailed in Turkheimer’s 2003 IEEE paper. This model requires input of the times, reftac and roitac and k2prime. Weights are optional.

First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.

taclow = tacdata$FSLSFC    # Frontal Cortex 
tacmed = tacdata$fslLST    # Limbic Striatum
tachigh = tacdata$fslAST   # Associative Striatum

refmlLogan_tstar(times, reftac, taclow, tacmed, tachigh, k2prime, 'demonstration')

Using this information (unshown), we see that we should probably use 9 included frames for the fit. We shall use 6.

refmlloganout <- refmlLogan(times, reftac, roitac, k2prime,
                        tstarIncludedFrames = 9, weights=weights)

Parameters

pandoc.table(refmlloganout$par, digits = 3)
bp k2
1.53 0.0437

Plot

The output object can be used to visualise the fit:

plot_refmlLoganfit(refmlloganout)
RefmlLogan Fit Check

RefmlLogan Fit Check

… and the residuals:

plot_residuals(refmlloganout)
RefLogan Residuals

RefLogan Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(refmlloganout))
  Length Class Mode
par 2 data.frame list
fit 13 lm list
tacs 3 data.frame list
fitvals 4 data.frame list
weights 14 -none- numeric
k2prime 1 -none- numeric
tstarIncludedFrames 1 -none- numeric
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(refmlloganout$fit), digits = 3)
Table continues below
r.squared adj.r.squared sigma statistic p.value df logLik AIC
1 1 94.4 48778 3.13e-15 2 -53.2 112
BIC deviance df.residual
113 62431 7

Patlak Reference Tissue Model

refPatlak requires input of the times, reftac and roitac. Weights are optional.

First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.

taclow = tacdata$FSLSFC    # Frontal Cortex 
tacmed = tacdata$fslLST    # Limbic Striatum
tachigh = tacdata$fslAST   # Associative Striatum
refPatlak_tstar(times, reftac, taclow, tacmed, tachigh, 'demonstration')

Since this is actually a reversible tracer, we shouldn’t really be using refPatlak at all, but we can use four frames to test it out (unshown).

refpatlakout <- refPatlak(times, reftac, roitac, 
                        tstarIncludedFrames = 4, weights=weights)

Parameters

pandoc.table(refpatlakout$par, digits=4)
## 
## --------
##    K    
## --------
## 0.005227
## --------

Plot

The output object can be used to visualise the fit:

plot_refPatlakfit(refpatlakout)
RefPatlak Fit Check

RefPatlak Fit Check

… and the residuals:

plot_residuals(refpatlakout)
RefPatlak Residuals

RefPatlak Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(refpatlakout))
  Length Class Mode
par 1 data.frame list
fit 13 lm list
tacs 3 data.frame list
fitvals 2 data.frame list
weights 14 -none- numeric
tstarIncludedFrames 1 -none- numeric
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(refpatlakout$fit), digits = 3)
Table continues below
r.squared adj.r.squared sigma statistic p.value df logLik AIC
0.894 0.84 0.0435 16.8 0.0547 2 8.06 -10.1
BIC deviance df.residual
-12 0.00378 2

Models Requiring Arterial Input

Data

First all the data is loaded (unshown). This data includes:

The data is plotted below:

par(mfrow=c(2,2))

plot(t_tac, tac)
plot(blooddata$Time.sec., blooddata$Cbl.nCi.cc., xlab = 'Time (sec)', ylab = 'Cb')
plot(plasmadata$Time.sec., plasmadata$Cpl.nCi.cc., xlab = 'Time (sec)', ylab = 'Cp')
plot(parentdata$Times, parentdata$Fraction, xlab = 'Time (sec)', ylab = 'Parent Fraction')
Compartmental Model Data

Compartmental Model Data

First, we must organise the blood data and put it together. This function can also accept parent-fraction-corrected plasma data too.

input <- blood_interp(t_blood = blooddata$Time.sec./60, 
                      blood = blooddata$Cbl.nCi.cc.,
                      t_plasma = plasmadata$Time.sec./60, 
                      plasma = plasmadata$Cpl.nCi.cc.,
                      t_parentfrac = parentdata$Times/60,
                      parentfrac = parentdata$Fraction)

One-Tissue Compartment Model

1TCM requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function and vB are also optional. If they are not included, then the function will fit them.

onetcmout <- onetcm(t_tac, tac, input, weights=weights)

Parameters

pandoc.table(onetcmout$par, digits = 3)
K1 k2 inpshift vB Vt
0.112 0.0317 -0.164 0.0824 3.54

Input Delay

The inpshift is the fitted delay. We can visualise the matching of the two curves without any time shifting (i.e. as they are in the same time scale):

plot_inptac_timings(t_tac, tac, input, inpshift = 0)
Delay Check Before

Delay Check Before

The timing of these two clearly do not match. Now we can visualise the fit with the time shifting fitted by the one tissue compartment model: in this case, the delay is negative, which means that the TAC was shifted to be a little bit later.

plot_inptac_timings(t_tac, tac, input, inpshift = onetcmout$par$inpshift)
Delay Check After

Delay Check After

This parameter has clearly performed a good job in fitting the timing of the two curves to one another.

Plot

The output object can be used to visualise the model fit:

plot_1tcmfit(onetcmout)
1TCM Fit Check

1TCM Fit Check

… and the residuals:

plot_residuals(onetcmout)
1TCM Residuals

1TCM Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(onetcmout))
  Length Class Mode
par 5 data.frame list
par.se 5 data.frame list
fit 6 nls list
tacs 3 data.frame list
input 4 data.frame list
weights 38 -none- numeric
inpshift_fitted 1 -none- logical
vB_fitted 1 -none- logical
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(onetcmout$fit), digits=3)
sigma isConv finTol logLik AIC BIC deviance df.residual
22.2 TRUE 1.49e-08 -182 373 381 14759 30

Two-Tissue Compartment Model

2TCM requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function and vB are also optional. If they are not included, then the function will fit them. In this case, we can use the delay from the 1tcm for the 2tcm fit, and it will not be fitted.

twotcmout <- twotcm(t_tac, tac, input, weights=weights, inpshift = onetcmout$par$inpshift)

Parameters

pandoc.table(twotcmout$par, digits = 3)
K1 k2 k3 k4 vB Vt
0.128 0.0721 0.0436 0.0295 0.0503 4.42

Plot

The output object can be used to visualise the model fit:

plot_2tcmfit(twotcmout)
2TCM Fit Check

2TCM Fit Check

… and the residuals:

plot_residuals(twotcmout)
2TCM Residuals

2TCM Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(twotcmout))
  Length Class Mode
par 6 data.frame list
par.se 6 data.frame list
fit 6 nls list
tacs 3 data.frame list
input 4 data.frame list
weights 38 -none- numeric
inpshift_fitted 1 -none- logical
vB_fitted 1 -none- logical
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(twotcmout$fit), digits=3)
sigma isConv finTol logLik AIC BIC deviance df.residual
7.73 TRUE 1.49e-08 -141 294 303 1731 29

Logan Plot

The Logan Plot method requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function and vB are also optional. If they are not included, then they will be set to zero. For the input delay, we can use the delay from the 1tcm fit. If vB is included, it will be corrected for before the parameter estimation using the following formula:

\[ C_{T}(t) = \frac{C_{Measured}(t) - vB\times C_{B}(t)}{1-vB} \]

First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.

lowroi = inptacdata$PUT
medroi = inptacdata$GM
highroi = inptacdata$gmCER

Logan_tstar(t_tac, lowroi, medroi, highroi, input, filename='demonstration', 
            inpshift = onetcmout$par$inpshift, vB = 0.05, gridbreaks=4)

Using this information (unshown), we see that we should probably use 10 frames for the fit.

loganout <- Loganplot(t_tac, tac, input, 10, weights, inpshift = onetcmout$par$inpshift, vB=0.05)

Parameters

loganout$par
##         Vt
## 1 4.225652

Plot

The output object can be used to visualise the model fit:

plot_Loganfit(loganout)
Logan Plot Fit Check

Logan Plot Fit Check

… and the residuals:

plot_residuals(loganout)
Logan Plot Residuals

Logan Plot Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(loganout))
  Length Class Mode
par 1 data.frame list
par.se 1 data.frame list
fit 13 lm list
tacs 2 data.frame list
fitvals 2 data.frame list
input 4 data.frame list
weights 38 -none- numeric
inpshift 1 -none- numeric
vB 1 -none- numeric
tstarIncludedFrames 1 -none- numeric
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(loganout$fit), digits=3)
Table continues below
r.squared adj.r.squared sigma statistic p.value df logLik AIC
0.998 0.998 1.49 4754 2.18e-12 2 -17.2 40.3
BIC deviance df.residual
41.2 17.8 8

MA1

The MA1 method requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function is also optional. If it is not included, then it will be set to zero.

First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.

lowroi = inptacdata$PUT
medroi = inptacdata$GM
highroi = inptacdata$gmCER

ma1_tstar(t_tac, lowroi, medroi, highroi, input, filename='demonstration', 
          inpshift = onetcmout$par$inpshift, gridbreaks=4)
## png 
##   2

MA1 tstar finder

Using this information, we see that we should probably use very few frames for the fit: 11 seems reasonably good and reasonably stable.

ma1out <- ma1(t_tac, tac, input, 11, weights, inpshift = onetcmout$par$inpshift)

Parameters

ma1out$par
##         Vt
## 1 4.143303

Plot

The output object can be used to visualise the model fit:

plot_ma1fit(ma1out)
MA1 Plot Fit Check

MA1 Plot Fit Check

… and the residuals:

plot_residuals(ma1out)
MA1 Plot Residuals

MA1 Plot Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(ma1out))
  Length Class Mode
par 1 data.frame list
fit 13 lm list
tacs 2 data.frame list
fitvals 6 data.frame list
input 4 data.frame list
weights 38 -none- numeric
inpshift 1 -none- numeric
tstarIncludedFrames 1 -none- numeric
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(ma1out$fit), digits=3)
Table continues below
r.squared adj.r.squared sigma statistic p.value df logLik AIC
0.999 0.999 6.6 4202 4.28e-14 2 -35.4 76.7
BIC deviance df.residual
77.9 392 9

MA2

The MA2 method requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function is also optional. If it is not included, then it will be set to zero.

ma2out <- ma2(t_tac, tac, input, weights, inpshift = onetcmout$par$inpshift)

Parameters

pandoc.table(ma2out$par, digits = 3)
Vt Vs Vnd K1 k2 k3 k4
4.11 3.42 0.684 0.153 0.223 0.00187 0.000373

One can also calculate the rate constants from the fit terms. The article provides a formula for calculation of VS, but they left out a number though for which term to use (this error is also on the PMOD page, which simply displays an image of the formula from the article). I believe I have chosen the correct one, but this might not be the case. VND can be calculated by subtracting Vs from VT.

Plot

The output object can be used to visualise the model fit:

plot_ma2fit(ma2out)
MA2 Plot Fit Check

MA2 Plot Fit Check

… and the residuals:

plot_residuals(ma2out)
MA2 Plot Residuals

MA2 Plot Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(ma2out))
  Length Class Mode
par 7 data.frame list
fit 13 lm list
tacs 2 data.frame list
fitvals 8 data.frame list
input 4 data.frame list
weights 38 -none- numeric
inpshift 1 -none- numeric
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(ma2out$fit), digits=3)
Table continues below
r.squared adj.r.squared sigma statistic p.value df logLik AIC
0.977 0.974 40.6 319 4.09e-24 4 -172 355
BIC deviance df.residual
363 49558 30

Multilinear Logan Plot

\[\int _{0} ^{t} C_{T} (\tau) d\tau = V_{T}\int _{0} ^{t} C_{P} (\tau) d\tau-\frac{1}{K_{2}}C_{T}(t)\]

The multilinear Logan Plot is the kinetic model used in WAPI, which is detailed in Turkheimer’s 2003 IEEE paper. This model requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function is also optional, and if it is not included, it will be set to zero.

First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.

lowroi = inptacdata$PUT
medroi = inptacdata$GM
highroi = inptacdata$gmCER

mlLogan_tstar(t_tac, lowroi, medroi, highroi, input, filename='demonstration', 
            inpshift = onetcmout$par$inpshift, gridbreaks=4)

Using this information (unshown), we see that we should probably use 10 frames for the fit.

mlloganout <- mlLoganplot(t_tac, tac, input, 10, weights, inpshift = onetcmout$par$inpshift)

Parameters

pandoc.table(mlloganout$par, digits=3)
Vt k2
4.21 0.0179

Plot

The output object can be used to visualise the model fit:

plot_mlLoganfit(mlloganout)
mlLogan Plot Fit Check

mlLogan Plot Fit Check

… and the residuals:

plot_residuals(mlloganout)
mlLogan Plot Residuals

mlLogan Plot Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(mlloganout))
  Length Class Mode
par 2 data.frame list
fit 13 lm list
tacs 2 data.frame list
fitvals 4 data.frame list
input 4 data.frame list
weights 38 -none- numeric
inpshift 1 -none- numeric
tstarIncludedFrames 1 -none- numeric
model 1 -none- character

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(mlloganout$fit), digits=3)
Table continues below
r.squared adj.r.squared sigma statistic p.value df logLik AIC
1 1 278 14762 5.39e-15 2 -69.4 145
BIC deviance df.residual
146 617196 8

Patlak Plot

The Patlak Plot method requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function is also optional. If it is not included, then it will be set to zero. For the input delay, we can use the delay from the 1tcm fit.

First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.

Patlak_tstar(t_tac, lowroi, medroi, highroi, input, filename='demonstration', 
            inpshift = onetcmout$par$inpshift, gridbreaks=4)

Since this is actually a reversible tracer, we shouldn’t really be using the Patlak Plot at all, but we can use nine frames to test it out (unshown).

patlakout <- Patlakplot(t_tac, tac, input, 9, weights, inpshift = onetcmout$par$inpshift)

Parameters

pandoc.table(patlakout$par, digits=3)
K
0.0118

Plot

The output object can be used to visualise the model fit:

plot_Patlakfit(patlakout)
Patlak Plot Fit Check

Patlak Plot Fit Check

… and the residuals:

plot_residuals(patlakout)
Patlak Plot Residuals

Patlak Plot Residuals

Output structure:

Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(patlakout))
  Length Class Mode
par 1 data.frame list
par.se 1 data.frame list
fit 13 lm list
tacs 2 data.frame list
fitvals 2 data.frame list
input 4 data.frame list
weights 38 -none- numeric
inpshift 1 -none- numeric
tstarIncludedFrames 1 -none- numeric

This means that the fit can be easily evaluated and compared with other fits:

pandoc.table(glance(patlakout$fit), digits=3)
Table continues below
r.squared adj.r.squared sigma statistic p.value df logLik AIC
0.997 0.997 0.0828 2772 2.34e-10 2 10.7 -15.4
BIC deviance df.residual
-14.8 0.048 7

Other Models

SIME

The SIME method requires input of the TAC times, several TACs, input and a grid of VND values. Weights are optional. ROI sizes can optionally be included to weight the contributions of the different ROIs. The time shift of the input function and vB are also optional. If they are not included, they will be set to 0 and 0.05 respectively.

tacdf <- data.frame(GM = inptacdata$GM, 
                    ACC = inptacdata$gmACC, 
                    CER = inptacdata$gmCER, 
                    PUT = inptacdata$PUT)

Vndgrid <- seq(from=0, to=5, by=0.1)

SIMEout <-  SIME(t_tac, tacdf, input, Vndgrid, weights = weights, 
    inpshift = onetcmout$par$inpshift, vB = 0.05, twotcmstart = 1)

Parameters

pandoc.table(SIMEout$par, digits=3)
Vnd
1.8

Plot

The output object can be used to visualise the model fit:

plot_SIMEfit(SIMEout)
SIME Plot Fit Check

SIME Plot Fit Check

Output structure:

Included in the output are all the parameters, the parameters used for the fitting functions, the fits themselves and the raw data.

pandoc.table(summary(SIMEout))
  Length Class Mode
par 1 data.frame list
tacs 7 data.frame list
fitvals 2 data.frame list
roifits 5 tbl_df list
input 4 data.frame list
weights 38 -none- numeric
roiweights 4 -none- numeric
inpshift 1 -none- numeric
vB 1 -none- numeric
model 1 -none- character

Alternatively, this model can be implemented in a faster manner of sequentially refining its results:

Vndgrid <- seq(from=0, to=5, by=1)

SIMEout1 <-  SIME(t_tac, tacdf, input, Vndgrid, weights = weights, 
    inpshift = onetcmout$par$inpshift)

Vndgrid <- seq(from=SIMEout1$par$Vnd-1, to=SIMEout1$par$Vnd+1, length.out = 10)

SIMEout2 <-  SIME(t_tac, tacdf, input, Vndgrid, weights = weights, 
    inpshift = onetcmout$par$inpshift)

Vndgrid <- seq(from=SIMEout2$par$Vnd-0.2, to=SIMEout2$par$Vnd+0.2, length.out = 10)

SIMEout3 <-  SIME(t_tac, tacdf, input, Vndgrid, weights = weights, 
    inpshift = onetcmout$par$inpshift)